Case Study 2 - Arabidopsis Time Series
Luke Yates
2025-02-18
About
This tutorial demonstrates the R package cpam for the
analysis of time series omics data. It reproduces the results for the
second case study presented in the manuscript Shape-constrained,
changepoint additive models for time series omics data with
cpam.
Data
The data are oublically available
Installation
You can install the package from GitHub using:
Getting started
Load packages
Experimental design and transcript-to-gene mapping
First we create the experimental design which must have at least the following columns: time, sample, and path. The rep column is optional and is used here to generate the sample names. We have already run kallisto to quantify transcript abundances (with 100 bootstraps), and here the path column contains the path to the abundance file for each sample.
ed <-
expand_grid(time = c(0,30,60,67.5,75,90,120), rep = 1:3) %>%
mutate(sample = paste0("pc_t",time,"_r",rep),
path = paste0("case_studies/crisp/data/kallisto/",sample,"/abundance.h5"))
edThe transcript-to-gene mapping for Arabidopsis thaliana can be downloaded from TAIR (…). We have downloaded it already and we load it now
This file should have two columns, target_id and gene_id.
Fitting the models
To fit the models, we first prepare the cpam object, then compute p-values, estimate the changepoint, and select the shape for each transcript. The last step takes the longest (here just under 13 minutes) but it is worth the wait to be able to visualise and cluster the transcripts by shape.
cpo <- prepare_cpam(exp_design = ed,
model_type = "case-only",
t2g = t2g,
import_type = "kallisto",
num_cores = 4)
cpo <- compute_p_values(cpo) # 1:52
cpo <- estimate_changepoint(cpo) # 6:32 secs
cpo <- select_shape(cpo) # 12:54 secs
We can look at a summary of the fitted cpam object
cpo
## cpam object
## -----------
## case-only time series
## 21 samples
## 7 time points
## Overdispersion estimated using 100 inferential replicates
## Counts rescaled by estimated overdispersionAlthough it cannot be launched here on the online tutorial, you can
launch the Shiny app on your own computer to visualise the results
interactively using visualise(cpo).
Result tables
The results of the analysis are summarised using the
results function.
The generated results can be filtered by specifying minimum counts, minimum log-fold changes, and maximum p-values. For example, to return only the transcripts with a log-fold change greater than 2, at least 50 counts, and a p-value less than 0.01, we can run
Plotting genes and transcripts
A single gene can be plotted using the plot_cpam
function. Here we plot the gene AT1G64140
The
subtitle shows (0,tp) indicating a changepoint at time point 0 (i.e., no
changepoint) and an unconstained ‘tp’ (thinplate) shape. This selection
of ‘tp’ suggests that the trend for this gene does not conform to one of
the simpler shape types that cpam uses. To force the
cpam to choose among the simpler forms, we set
shape_type = "shape2" in the plot_cpam
function. For example:
Here
a concave shape (‘cv’) is chosen, and we can see this trend deviates
from the data substantially more that the unconstrained shape.
Next we plot a gene with multiple transcripts, some of which have changepoints different from 0.
The first transcript has a changepoint at 30 mins and the second at 0.
Both have an unconstrained shape. The transcripts can plotted separately
by setting facet =T in the plot_cpam function.
There are may other settings that can be adjusted in the
plot_cpam function, see the documentation for more details.
Changepoints, shapes and other results from the fitted models can also
be extracted manually from the cpam object. For example, to
extract the shape of the transcripts
Lastly, we plot the two remaining genes shown in the manuscript, AT4G34590 and AT3G23280.
Clusters
The results function can be used to generate clusters according to
selected filters. Here we generate clusters of transcripts with at least
100 counts, a log-fold change greater than 1, and a p-value less than
0.01. The plot_cluster function can then be used to
visualise the clusters which we do here for targets with changepoints at
30 mins and the ‘micv’ (montonic increasing concave) shape.
res <- results(cpo, min_count = 100, min_lfc = 1, p_threshold = 0.01)
plot_cluster(cpo, res, changepoints = 30, shapes = "micv")
Clustering can be further refined based on, for example, the rate at
which the above transcripts attain their maximum values. We illustrate
advanced refinements such as this in our other case study here.
Advanced plotting
OSE rule plots
When estimating the changepoint, cpam applies the
one-standard error to mitigate overfitting by taking into account model
selection uncertainty. To visualise the rule for a given target, we plot
the pointwise score differences and their standard errors for all
changepoints. First we extract the score table which contains the
pointwise cross-validation scores for each model (changepoint).
score_table <-
cpo$changepoints %>%
filter(target_id == "AT3G23280.1") %>%
pull(score_table) %>%
.[[1]]
score_tableNext, we determine the minimum scoring model, and compute the pointwise score differences.
m.min <- score_table %>% purrr::map_dbl(mean) %>% which.min() %>% names; m.min
## [1] "30"
score_table_diff <- score_table - score_table[[m.min]]
score_table_diffFinally, we compute and plot the mean score differences and their standard errors. The OSE-selected changepoint is latest changepoint whose score is within one standard error of the minimum score difference.
tibble(cp = names(score_table_diff) %>% {factor(.,.)},
score_diff = score_table_diff %>% purrr::map_dbl(mean),
se_diff = score_table_diff %>% purrr::map_dbl(~sd(.x)/sqrt(nrow(score_table_diff)))) %>%
ggplot(aes(x = cp)) +
geom_linerange(aes(ymin = score_diff - se_diff,
ymax = score_diff + se_diff)) +
geom_hline(yintercept = 0, lty = "dashed") +
geom_point(aes(y = score_diff, colour = "min"), shape = 1, size = 5,
data = ~ filter(.x, score_diff == min(score_diff))) +
geom_point(aes(y = score_diff, colour = "1se"), shape = 1, size = 5,
data = ~ filter(.x, se_diff >= score_diff) %>%
filter(as.numeric(cp) == max(as.numeric(cp)))) +
geom_point(aes(y = score_diff), size = 2) +
scale_colour_manual(values = c(`1se` = "#CB181D", min = "#08519C"),
breaks = c("min","1se")) +
labs(col = NULL, y = expression(Delta*"P"), x = "Changepoint",
subtitle = "One-standard-error rule") +
theme_classic() +
theme(legend.position = "bottom") 
For
the target AT3G23280.1, the OSE rule selects the changepoint at 67.5
mins, although we can see that the subsequent point at 75 mins is a
close contender (compare the OSE plot with the plot of data and fitted
model). For comparison, we plot the fitted model using the minimum
scoring model instead of the OSE-selected model.
plot_cpam(cpo, target_id = "AT3G23280.1", cp_type = "cp_min", bs = "tp") +
labs(subtitle = "Minimum scoring model")
In
general, the minimum scoring model captures a more complex trend than
the OSE-selected model, although the OSE-selected model is likely to
generalise better to new data. See the manuscript for a simulation study
comparing two selection rules.
Isoform structure plots
Here we show how to reproduce the isoform structure plots for the gene AT3G23280.
First, download the the gtf file for Arabidopsis thaliana from
Ensembl Plants. Using the rtracklayer package, import the
gtf file and convert it to a tibble.
# Download the file
download.file(
url = "https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-55/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.55.gtf.gz",
destfile = "Arabidopsis_thaliana.TAIR10.55.gtf.gz",
method = "auto"
)
gtf <- rtracklayer::import("Arabidopsis_thaliana.TAIR10.55.gtf.gz") %>% dplyr::as_tibble()This next part is a bit more involved and requires some manual adjustments to achieve the desired result.
library(ggtranscript)
g_annotation <- gtf %>% filter(gene_id == "AT3G23280")
g_annotation_exons <- g_annotation %>% dplyr::filter(type == "exon")
g_annotation_cds <- g_annotation %>% dplyr::filter(type == "CDS")
tx_ids <- unique(g_annotation_exons$transcript_id)
g_annotation_exons_1 <- g_annotation_exons %>% filter(transcript_id == tx_ids[1])
g_annotation_exons_2 <- g_annotation_exons %>% filter(transcript_id == tx_ids[2])
g_annotation_cds_1 <- g_annotation_cds %>% filter(transcript_id == tx_ids[1])
g_annotation_cds_2 <- g_annotation_cds %>% filter(transcript_id == tx_ids[2])
tx_cols <- c("#CB181D","#08519C")
p_exon <-
g_annotation_cds_1 %>%
ggplot(aes(
xstart = start,
xend = end,
y = "y"
)) +
geom_half_range(
aes(fill = transcript_id, col = transcript_id),
height = 0.125,
range.orientation = "top",
alpha = 0.1,
linewidth = 0.5,
position = position_nudge(y = 0.001)
) +
geom_half_range(
aes(fill = transcript_id, col = transcript_id),
data = g_annotation_cds_2,
height = 0.125,
alpha = 0.1,
linewidth = 0.5,
position = position_nudge(y = -0.001)
) +
geom_intron(
data = to_intron(g_annotation_exons_1, "transcript_id"),
col = "grey20",
arrow = grid::arrow(ends = "last", length = grid::unit(2, "mm")),
linewidth = 0.5,
arrow.min.intron.length = 150,
strand = g_annotation_cds_1$strand[1]
) +
scale_fill_manual(values = tx_cols, aesthetics = c("color","fill")) +
theme_void() +
theme(
plot.margin = unit(c(-1.5, 1, -1.5, 1), "cm"),
legend.position = "none"
) +
annotate("segment", x = min(g_annotation_cds_2$start), xend = max(g_annotation_cds_2$end),
y = "y", yend = "y", col = "grey20", linewidth = 0.8) +
xlim(min(g_annotation_cds_2$start) - 1200,NA) +
annotate("text", x = min(g_annotation_cds_2$start) - 50, y = 1 + 0.125/2, label = tx_ids[1],
hjust = 1, size = 3.5, col = tx_cols[1]) +
annotate("text", x = min(g_annotation_cds_2$start) - 50, y = 1 - 0.125/2, label = tx_ids[2],
hjust = 1, size = 3.5, col = tx_cols[2]);p_exon
Extract plot data for full customisation
If you require more customisation than the plot_cpam
function allows, you can extract the model fits and plot their
predictions manually. Here we show how to extract the model fits, model
predictions, and observed data for a given target.
# extract the model fit(s)
fit <- plot_cpam(cpo, target_id = "AT3G23280.1",return_fits_only = T)
# generate model prediction data
cpam:::predict_cpgam(fit, logged = F, length.out = 200, ci_prob = "se")The count data can be extract the cpam object
Session Info
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Pop!_OS 22.04 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-openmp/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Australia/Hobart
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggtranscript_1.0.0 ggplot2_3.5.1 stringr_1.5.1 tidyr_1.3.1 readr_2.1.5 dplyr_1.1.4
## [7] cpam_0.0.0.9000
##
## loaded via a namespace (and not attached):
## [1] SummarizedExperiment_1.36.0 gtable_0.3.6 rjson_0.2.23 xfun_0.50
## [5] bslib_0.9.0 ggrepel_0.9.6 Biobase_2.66.0 lattice_0.22-6
## [9] tzdb_0.4.0 vctrs_0.6.5 tools_4.4.1 bitops_1.0-9
## [13] generics_0.1.3 parallel_4.4.1 stats4_4.4.1 curl_6.2.0
## [17] tibble_3.2.1 pkgconfig_2.0.3 Matrix_1.7-2 RColorBrewer_1.1-3
## [21] S4Vectors_0.44.0 lifecycle_1.0.4 GenomeInfoDbData_1.2.13 scam_1.2-18
## [25] compiler_4.4.1 farver_2.1.2 Rsamtools_2.22.0 Biostrings_2.74.1
## [29] munsell_0.5.1 codetools_0.2-20 GenomeInfoDb_1.42.3 htmltools_0.5.8.1
## [33] sass_0.4.9 RCurl_1.98-1.16 yaml_2.3.10 crayon_1.5.3
## [37] pillar_1.10.1 jquerylib_0.1.4 BiocParallel_1.40.0 DelayedArray_0.32.0
## [41] cachem_1.1.0 abind_1.4-8 nlme_3.1-167 aggregation_1.0.1
## [45] tidyselect_1.2.1 digest_0.6.37 stringi_1.8.4 purrr_1.0.4
## [49] restfulr_0.0.15 labeling_0.4.3 splines_4.4.1 fastmap_1.2.0
## [53] grid_4.4.1 SparseArray_1.6.1 colorspace_2.1-1 cli_3.6.4
## [57] magrittr_2.0.3 S4Arrays_1.6.0 XML_3.99-0.18 withr_3.0.2
## [61] scales_1.3.0 UCSC.utils_1.2.0 rmarkdown_2.29 XVector_0.46.0
## [65] httr_1.4.7 matrixStats_1.5.0 hms_1.1.3 evaluate_1.0.3
## [69] knitr_1.49 GenomicRanges_1.58.0 IRanges_2.40.1 BiocIO_1.16.0
## [73] rtracklayer_1.66.0 mgcv_1.9-1 rlang_1.1.5 Rcpp_1.0.14
## [77] glue_1.8.0 BiocGenerics_0.52.0 rstudioapi_0.17.1 jsonlite_1.8.9
## [81] R6_2.6.1 MatrixGenerics_1.18.1 GenomicAlignments_1.42.0 zlibbioc_1.52.0